Thermodynamically consistent mesoscopic fluid particle models for a van der Waals 
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The GENERIC structure allows for a unified treatment of different discrete models of hydrody- 
namics. We first propose a finite volume Lagrangian discretization of the continuum equations of 
hydrodynamics through the Voronoi tessellation. We then show that a slight modification of these 
discrete equations has the GENERIC structure. The GENERIC structure ensures thermodynamic 
consistency and allows for the introduction of correct thermal noise. In this way, we obtain a con- 
sistent discrete model for Lagrangian fiuctuating hydrodynamics. For completeness, we also present 
the GENERIC versions of the Smoothed Particle Dynamics model and of the Dissipative Particle 
Dynamics model. The thermodynamic consistency endorsed by the GENERIC framework allows 
for a coherent discussion of the gas-liquid phase coexistence of a van der Waals fiuid. 



I. INTRODUCTION 

The behaviour of complex fluids like colloids, emul- 
sions, polymers or multiphasic fluids is affected by the 
strong coupling between the microstructure of these 
fluids with the macroscopic flow. The complexity of 
these systems requires the use of novel computer sim- 
ulations techniques and algorithms. Usual macroscopic 
approaches that solve partial differential equations with 
constitutive equations are not very useful because the ba- 
sic input, the constitutive equation is usually not known. 
Also, these approaches neglect the presence of thermal 
noise, which is the responsible for the Brownian motion 
of small suspended objects and, therefore, for the diffu- 
sive processes that affect the microstructure of the fluid. 
In recent years, there has been a large effort in order 
to develop mesoscopic techniques in order to tackle the 
problems arising in the simulation of complex fluids. 

Dissipative Particle Dynamics is a mesoscopic parti- 
cle based simulation method that allows one to model 
hydrodynamic behavior and that it consistently includes 
thermal fluctuations. It was introduced by Hoogerbrugge 
and Koelman in 1991 under the motivation of designing 
an off-lattice algorithm inspired by the ideas behind the 
lattice-gas method [Q. Since then, the model has re- 
ceived a great deal of attention. From a theoretical point 
of view, the model has been given a solid background as 
a statistical mechanics model . The hydrodynamic be- 
havior has been analyzed and the methods of kinetic 
theory have provided explicit formulae for the transport 
coefficients in terms of the model parameters ^ . A gen- 
eralization of DPD has also been presented in order to 
conserve energy From the side of applications, the 
method is very versatile and has proven to be useful in 
the simulation of flows in porous media Q , colloidal sus- 
pensions 1^,01, polymer suspensions microphase sep- 
aration of copolymers , multicomponent flows |ic| ] and 
thin-fllm evolution lllll. 



The physical picture behind the dissipative particles 
used in the model is that they represent mesoscopic por- 
tions of real fluid, say clusters of molecules moving in 
a coherent and hydrodynamic fashion. The interaction 
between these particles are postulated from simplicity 
and symmetry principles that ensure the correct hydro- 
dynamic behavior. DPD faces, however, a conceptual 
problem. The thermodynamic behavior of the model 
is determined by the conservative forces introduced in 
the model. This forces are assumed to be soft forces in 
counter distinction to the singular forces of the Lennard- 
Jones type used in molecular dynamics. But there is 
no well-deflned procedure to relate the shape and ampli- 
tude of the conservative forces with a prescribed thermo- 
dynamic behavior (although attempts in that direction 
have been undertaken, see Ref. |^). Also, it is not clear 
which physical time and length scales the model actually 
describes, even though the presence of thermal noise sug- 
gests the foggy area of the mesoscopic realm. We will see 
that both problems are closely related. 

Dissipative Particle Dynamics is very similar in spirit 
to the popular method of Smoothed Particle Dynamics 
The method was introduced in the context of as- 
trophysics computation in the early 70's [|l2| and very 
recently it has been applied to the study of laboratory 
scale viscous and thermal flows |14| in simple geome- 
tries. SPD is essentially a Lagrangian discretization of 
the Navier-Stokes equations by means of a weight func- 
tion. The procedure transforms the partial differential 
equations of continuum hydrodynamics into ordinary dif- 
ferential equations. These equations can be further in- 
terpreted as the equations of motion for a set of particles 
interacting with prescribed laws of force. The technique 
thus allows one to solve PDE's with molecular dynam- 
ics codes. Again, these particles can be understood as 
physical portions of the fluid that evolve coherently along 
the flow. The problem with Smoothed Particle Hydro- 
dynamics is that it does not include thermal fluctuations 
and, therefore, cannot be applied to the study of complex 
fluids at mesoscopic scales. 
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We have recently shown that the conceptual problems 
in DPD and the inclusion of thermal fluctuations in SPD 
can be resolved by formulating convenient generaliza- 
tions of both methods under the general framework of 
GENERIC ||T^. In the present paper, we take a further 
look at the problem of formulating consistent models for 
the simulation of hydrodynamic problems. Our point of 
view here is to construct a finite volume algorithm for 
solving the Navier-Stokes equations in such a way that 
thermodynamic consistency is retained: The resulting al- 
gorithm conserves mass, momentum, and energy, and the 
entropy is an increasing function of time. Most impor- 
tant, we show how to include thermal noise in a consistent 
way, this is, producing the Einstein distribution function. 
We end up, therefore, with an algorithm for simulat- 
ing fluctuating hydrodynamics in a Lagrangian way [|l6| . 
This algorithm can be used as the basis for simulating 
colloidal suspensions, where the Brownian motion of the 
colloidal particles is due to the thermal fluctuations on 
the solvent [0 . We show also another potential applica- 
tion to multiphasic flow of the gas-liquid type. When the 
fluid is described with a van der Waals equation of state, 
the Einstein equilibrium distribution allows to discuss the 
gas-liquid phase transition in probabilistic terms. Such 
probabilistic approaches to equilibrium phase transitions 
have been studied in the past iQ. We should note, how- 
ever, that the GENERIC finite volume algorithm should 
allow us to study fully non-equilibrium situations created 
by external boundary conditions that can drive the sys- 
tem out of equilibrium. Start up of boiling of water in 
a pot could be addressed with the proposed GENERIC 
finite volume algorithm. 

Other very recent approaches to the study of the flow 
of liquid-vapor coexisting fluids are the Lattice Boltz- 
mann model introduced by Swift, Osborn and Yeomans 
1^ and improved in Refs. (20) in order to have a thermo- 
dynamically consistent model for a dense fluid that may 
exhibit liquid-vapor coexistence. However, only isother- 
mal models have been considered up to now. A second 
very promising approach is the Direct Simulation Monte 
Carlo method that has been conveniently generalized 
to deal with dense liquids with liquid-vapor coexistence 

H- 

In order to construct the finite volume algorithm we 
have been strongly inspired by the work of Flekkoy and 
Coveney p3|| . In that paper, the authors present a 
"bottom-up" derivation of Dissipative Particle Dynam- 
ics. Physical space is divided into Voronoi cells and ex- 
plicit definitions for the mass, momentum and energy 
of the cells in terms of the microscopic degrees of free- 
dom (positions and momenta of the constituent molecules 
of the fiuid) are given. The time derivatives of these 
phase functions have the structure of "microscopic bal- 
ance equations" in a discrete form. These equations are 
then divided into "average" and "stochastic" parts. To 
further advance into the formulation of a practical algo- 



rithm, the authors then propose phenomenological, phys- 
ically sensible, expressions for the average part and re- 
quire the fulfillment of the fluctuation-dissipation theo- 
rem for the stochastic part. Because of the use of the 
phenomenological expressions, we cannot consider this a 
"bottom-up" approach. Strictly speaking, a bottom-up 
approach would require the use of a projection operator 
technique or kinetic theory, in order to relate the trans- 
port coefficients with the microscopic dynamics of the 
system (in the form of Grccn-Kubo formulae, for exam- 
pie). 

Instead, we propose in this paper a conspicuous "top- 
down" approach in which the deterministic continuum 
equations of hydrodynamics are the starting point. By 
making intensive use of the smooth Voronoi tessellation 
discovered by Flekkoy and Coveney the form of the dis- 
crete equations is dictated by the very structure of the 
continuum equations. Our approach is similar to that in 
Ref. [p4[ . However, we make a further requirement on 
the resulting finite volume discretization, which is that 
they must have the GENERIC structure. This enforces 
the addition of a tiny bit into the momentum equation. 
Having the GENERIC structure, it is trivial to obtain the 
stochastic part, which will be given by the fluctuation- 
dissipation theorem. In the concluding section we will 
show the similarities and differences between our equa- 
tions and those derived by Flekkoy and Coveney. 

The approach presented in this work has also a strong 
resemblance with Yuan and Doi simulation method which 
also uses the Voronoi tessellation in Lagrangian form |2^] 
(see also |26j). They have applied the method to the sim- 
ulation of concentrated emulsions under flow and several 
other applications to complex fluids are mentioned. The 
main difference between our work and that of Ref. p^ ] 
is the special care we have taken in order to have ther- 
modynamic consistency through the GENERIC formal- 
ism. This allows, among other things, to include correct 
thermal noise and describe diffusive aspects produced by 
Brownian motion on mesoscopic objects. Another differ- 
ence is that we deal with a compressible fluid in which the 
pressure is given through the equation of state as a func- 
tion of mass and entropy densities, in counterdistinction 
with Ref. where the pressure is obtained by satis- 
fying the incompressibility condition. Compressibility is 
necessary if gas-liquid coexistence is to be described. 
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II. REVIEW OF GENERIC 

For the sake of completeness we review in this sec- 
tion the GENERIC formahsm developed by Ottinger 
and Grmela |^ . The formalism of GENERIC (acronym 
for General Equation for Non Equilibrium Reversible Ir- 
reversible Coupling) states that all physically sounded 
transport equations in non-equilibrium thermodynamics 
have the same structure. A large body of evidence con- 
firms this assertion: Linear irreversible thermodynamics, 
non-relativistic and relativistic hydrodynamics, Boltz- 
mann's equation, polymer kinetic theory, and chemical 
reactions, just to mention a few, have all the GENERIC 
structure |£^,|9l. The GENERIC formalism is not only 
a way of rewriting known transport equations in a phys- 
ically transparent way, but it allows us to derive dynam- 
ical equations for new systems not considered so far in 
an astonishing simple way. The very structure of the for- 
malism ensures thermodynamic consistency, energy con- 
servation and positive entropy production. 

The essential assumption on which GENERIC is 
founded is that the relevant variables x used to describe 
the system at a certain level of description evolve in a 
time scale well-separated from the time scales of other 
variables in the system. In other words, the variables 
should provide a closed description in which the present 
state of the system depends only on the very recent past 
and memory effects can be neglected. This is a recurrent 
theme in non-equilibrium statistical mechanics since the 
pioneering works of Zwanzig and Mori . Actually, the 
GENERIC structure can be deduced from first principles 
by using standard projection operator formalism under a 
Markovian approximation p8| . 

Two basic building blocks in the GENERIC formalism 
are the energy E{x) and entropy S{x) functions of the 
variables x describing the state of the system at a par- 
ticular level of description H]. The GENERIC dynamic 
equations are given then by 



— - L — M — 

dt dx dx 



(1) 



The first term in the right hand side is named the re- 
versible part of the dynamics and the second term is 
named the irreversible part. The predictive power of 
GENERIC relies in the fact that very strong require- 
ments exists on the matrices L, M leaving small room 
for the physical input about the system. First, L is anti- 
symmetric whereas M is symmetric and positive semidef- 
inite. Most important, the following degeneracy condi- 
tions should hold 



L- 



dS_ 

dx 



0, 



M— = 0. 

ox 



(2) 



These properties ensure that the energy is a dynami- 
cal invariant, E = Q, and that the entropy is a non- 
decreasing function of time, 5 > 0, as can be proved by a 



simple application of the chain rule and the equations of 
motion (|l|). In the case that other dynamical invariants 
I{x) exist in the system (as, for example, linear or angu- 
lar momentum), then further conditions must be satisfied 
by L, M. In particular. 



— L—-0 

dx dx 



^ M —~0 
dx dx 



(3) 



which ensure that / = 0. 

The deterministic equations (0) are, actually, an ap- 
proximation in which thermal fluctuations are neglected. 
If thermal fluctuations are not neglected, the dynamics is 
described by stochastic differential equations or, equiv- 
alently, by a Fokker-Planck equation that governs the 
probability distribution function p = p{x,t). This FPE 
has the form l|2il 



dx 



L- 



dE 
dx 



dx 



dx 



(4) 



where ks is Boltzmann's constant. 

The distribution function of the variables of a given 
system at equilibrium is given by the Einstein distribu- 
tion function. This assertion can be proved under quite 
general hypothesis on the mixing character of the micro- 
scopic dynamics of the system [Q. If the microscopic 
dynamics ensures the existence of dynamical invariants 
like the energy E{x) and, perhaps, other invariants I{x), 
then the Einstein distribution function takes the form 



p'=^(x) = g{E{x),I{x)) exp{S{x)/kB}, 



(5) 



where the function g is completely determined by the 
arbitrary initial distribution of dynamical invariants. For 
example, if at an initial time the value of the invariants 
E{x),I(x) are known with high precision to be Eq,Iq, 
then the Einstein distribution function takes the form 



p^^ix) 



S{E{x) - Eo)S{I{x) ~ Iq) 
n{Eo,Io) 



e^kg'Six)}, (6) 



where f2(£'o:^o) is the normalization. Given the general 
argument behind the Einstein distribution function , 
it is sensible to demand that the Fokker-Planck equation 
(^ has as its (unique) equilibrium distribution function 
the Einstein distribution. This can be achieved, actu- 
ally, if the following further conditions on the form of the 
matrices L,M hold, 



d_ 

dx 



L- 



dE 
dx 



= 0, 



dx 



(7) 



The first property can be derived independently with pro- 
jection operator techniques psf whereas the second prop- 
erty implies that the last equation in (|^) is automatically 
satisfied. 
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When fluctuations are present, the entropy function 
S{x) might be a decreasing function of time. However, if 
one considers the entropy functional 



S{x)p{x,t)dx ~ ks / p{x, t) In p{x,t)dx, (8) 



S[pt 



it is possible to prove by using the Fokker-Planck equa- 
tion (Q) that dtS[pt] > 0. In other words, the entropy 
functional plays the role of a Lyapunov function. 

The stochastic differential equations that are mathe- 
matically equivalent to the above Fokker-Planck equation 
are given, with Ito interpretation, by [pll 



dx 



dE ^dS , (9 
L-— + M-— + kBj--M 
ox ox ox 



dt + dx. 



(9) 



to be compared with the deterministic equations . The 
stochastic term dx in Eqn. (^) is a linear combination of 
independent increments of the Wiener process. It satis- 
fies the mnemotechnical Ito rule 



dxdx"^ = 2kBMdt, 



(10) 



which means that dx is an infinitesimal of order 1/2 [ pi[ . 
Eqn. (|l^ is a compact and formal statement of the 
fluctuation-dissipation theorem. 

When formulating new models it might be convenient 
to specify dx directly instead of M. This ensures that M 
through (|lO|) automatically satisfies the symmetry and 
positive definite character. In order to guarantee that 
the total energy and dynamical invariants do not change 
in time, a strong requirement on the form of dx holds, 



——■dx — O, —-dx — O, 

Ox Ox 



(11) 



implying the last equations in (^) and (^. The geomet- 
rical meaning of ( pi] ) is clear. The random kicks pro- 
duced by dx on the state x are orthogonal to the gradi- 
ents of E, I. These gradients are perpendicular vectors 
(strictly speaking they are one forms) to the hypersurface 
E{x) = Eq, I{x) = /q. Therefore, the kicks let the state x 
always within the hypersurface of dynamical invariants. 

We finally close this section by noting that, formally, 
the size of the fluctuations is governed by the Boltzmann 
constant ks- If we take /cs ^ 0, the stochastic differen- 
tial equation (|^) becomes the deterministic equations (|l|) 
and the Fokker-Planck equation (^ has only first deriva- 
tives in the same way as the Liouville equation. In this 
limit, the distribution function p{x,t) does not show dis- 
persion and it is essentially a Dirac delta function eval- 
uated on the solution of the deterministic equation. In 
this case, the entropy functional S[p] reduces to the en- 
tropy S{x) and the entropy is a non-decreasing function 
of time. 



III. FINITE VOLUME METHOD WITH 
VORONOI CELLS 

As a first step for deriving the GENERIC equations 
for a model of fluid particles, we consider the method of 
finite volumes for the numerical integration of the equa- 
tions of continuum hydrodynamics. The basic motivation 
for this is to have a set of reference equations that serve 
as guidelines for the modeling in the GENERIC equa- 
tions in such a way that we can consider the GENERIC 
equations as reasonable approximations to the equations 
of continuum hydrodynamics. 

The finite volume method consists on integrating the 
continuum equations of hydrodynamics in a finite region 
of space (or finite volume) in such a way that ordinary 
differential equations for the average fields over the finite 
regions emerge. In this section we present a finite volume 
method that uses the Voronoi construction as a conceptu- 
ally and mathematically elegant method for discretizing 
the continuum equations of hydrodynamics. 

Following Flekkoy and Coveney [P3| , we introduce the 
smoothed characteristic function of the Voronoi cell /i 



X/.(i-) 



E.A(| 



(12) 



where the function A(r) = exp{— r^/2(T^} is a Gaussian 
of width a. When a ^ 0, the smoothed characteristic 
function tends to the actual characteristic function of the 
Voronoi cell, this is 



\im x^{v) =l[0i\v -K^l - \v -K, 



(13) 



where 9{x) is the Heaviside step function. The Voronoi 
characteristic function (^3|) takes the value 1 if r is nearer 
to than to any other H.i, with v ^ p.. Note that the 
characteristic function produces a covering of all space 
(i.e., a partition of unity), this is. 



(14) 



We introduce the volume of the Voronoi cell through 



JVt 



which satisfies the closure condition 



(15) 



(16) 



where Vt is the total volume. In Fig. |^ we show the 
Voronoi tessellation corresponding to 15 particles seeded 
at random in a periodic box. For an introduction to the 
Voronoi tessellation see WA- 
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We mention now some other useful properties of the 
smoothed characteristic function that will be needed 
later. First, due to the Gaussian form of A(r), 



VA(r) = — ^A(r)r 



Therefore, 
d 



dr 



rXM(i")(r 



By using the following property 



(17) 



(18) 



(19) 



which can be proved by using the definition (|T^) , one can 
rewrite Eqn. (18) as 



d 1 ^ 



(20) 



FIG. 1. Contour line at Xy-i^) = 0-5 for the set of all 
characteristic functions of the 15 particles located at random 
in a periodic box of size L = 100. The value of a is 0.03. 



Another useful relation is, 



d 1 

XmW = (5^^^XmW(i" - R-m) 



rXM(i")Xz^(i')(r-R^). 



(21) 



"Si 



We now introduce the following two quantities 
A^,^ = R^^ j ^Xm(i")X!^(i"), 



^LLU 



;Xt.ir)Xuir) r- 



(22) 



where i?^^ = |R^ — R^|. In appendix ^ we show that 
in the limit a the quantity A^^^ is actually the area 
of the contact face between Voronoi cells /i and v, 
whereas the vector c^i, is the position of the center of 
mass of the face fiv with respect to the "center" of the 
face (R^ +R^)/2. 



A. Balance equations 



FIG. 2. Contour lines of the function Xm('') f^i" particle 
= 1 of the previous tessellation. Inside the closed region 
the value of XM(r) is 1 and outside is zero. 



We introduce now the cell average [</>]^(t) over the 
Voronoi cell of an arbitrary density field 0(r, t) 



dr0(r,i)xp(r). 



(23) 



We will refer to [(?!']^(t) as a cell variable and we will see 
that it is an approximation for the value of the field at 
the discrete points given by the cell centers. 
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In principle, the Voronoi cell centers are allowed to 
move in an arbitrary way, this is, R^(t) are prescribed 
functions of time. The time derivative of the cell averages 
is given by 



+ 



1 



drx^{t)dt(i>ir,t), 



(24) 



where the dot means the time derivative. We see that 
[(j)]^{t) changes due to both, the motion of the cells and 
the intrinsic dependence of the field (j){r, t) on time. 

Now, we will assume that the field (/)(r, t) obeys a bal- 
ance equation, this is 



5t0(r,t) = -V-J(r,t), 



(25) 



where J(r,f) is an appropriate current density. By inte- 
gration by parts of the nabla operator and use of Eqns. 
( pO| ) and ( |l9| ) one arrives easily at the following expres- 
sion 



[4>Ut) 



(26) 



where 



R 



and we have introduced the face averages 



(27) 



ill ^ 



dr 



R,, + R, 



(28) 



Note that in the limit of sharp boundaries ti — > 0, is 
a vector which is parallel to the face fiv, whereas e^j^ is 
perpendicular to the face. 



We can write Eqn. (26) in the form 



■fiiy 



[J] 



Ri/ 



A 



(29) 



which satisfies 



(30) 



due to the symmetries [■ • -Jp,, = [• • -U,, [■ • -Jil^ = [■ • -jU^ 
of the face averages. Equation (pH) shows that the 
Voronoi discretization of the balance equation (|2^) con- 
serves exactly the extensive variables (which are of the 
form density X volume). 



B. Finite volumes for an inviscid fluid 

In what follows we will apply the method of finite vol- 
umes to the continuum equations of hydrodynamics. For 
the sake of clarity, we first consider the reversible part 
of the equations, which correspond the the usual Euler 
equations for an inviscid fluid. In the next subsection we 
consider the irreversible part. 

The Euler equations for an inviscid fluid are the con- 
tinuity equation 



dtp{v,t) - -V-g(r,i), 
the momentum balance equation 



(31) 



atg(r, t) = -VP(r, t) - V-g(r, i)v(r, t), (32) 
and the entropy equation 



dts{Y,t) = -V-s(r,i)v(r,t). 



(33) 



In these equations, p(r, t) is the mass density field, 
g(r,i) = p(r, i)v(r,t) is the momentum density field, 
with v(r, i) the velocity field, and s(r,t) is the entropy 
density field (entropy per unit volume). The pressure 
field P(r, t) is given, according to the local equilib- 
rium assumption, by P{v,t) = P°'i(p(r, i), s(r, i)) where 
P°'^{p, s) is the equilibrium equation of state that gives 
the macroscopic pressure in terms of the mass and en- 
tropy densities. 

We now write Eqn. (29) for the case that (j) = p. If the 



Voronoi cells do not move, then R^ = and expression 



tj simplifies considerably. 



dt 



(34) 



where A/^ ~ Vf^lp]^ is the mass of the Voronoi cell p. 
This corresponds to a Eulerian discretization of the con- 
tinuity equation (|3l|). The physical meaning of Eqn. ( |3^ ) 
is clear: A^^e^^-lpv]^^ is the total mass per unit time 
that crosses face piy and the rate of change of the mass 
of the cell p is the sum of this quantity for each face pv. 

The Lagrangian discretization of the continuity equa- 
tion in Eqn. ( |3l| ) is obtained by specifying the motion of 
the cells according to 
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(35) 



In this way, the Voronoi cells "follow" (the best they can) 
the flow field. We can write Eqn. (^9|) for the case ( |3l| ) 
as 



M„ 



A 



(36) 



The momentum balance equation (|3^) can similarly be 
treated and the following Lagrangian finite volume equa- 
tion arise 

V 



(37) 



where we have introduced the total momentum of cell /i 



Finally, the entropy equation Eqn. (|3^) takes the fol- 
lowing Lagrangian discretization 



where 5^ = V^s^ is the total entropy of cell 

In all these evolution equations, the subscript |rcv de- 
notes that the equations are actually the reversible part 
of the dynamics of a truly viscous fluid. 



C. Gradient expansion 

The above equations (^6|), (|37|), and ( ^8|) are rigorous 
and exact and do not depend on the typical size of the 
Voronoi cells. By taking sufficiently small cells, we can 
approximate the equations and transform them into a 
closed set of equations for the cell variables. In this way, 
a computationally feasible algorithm can be proposed for 
the updating of cell variables. 

Let us assume that the hydrodynamic fields have a 
typical length scale of variation \h and that the typical 
distance between cell centers is Ac- We introduce the 
resolution parameter as r = Xc/Xh and will assume that 
this parameter is very small. Therefore the dimensionless 
quantity 



A. 



(Rm) 



V0(R^) 



(39) 



will be very small. By denoting with C'(V) those terms 
which are of relative size r, we can write the cell average 
as 

1 



drx^(r)0(r) 



V, 



^ J (irx;.(r)0(r-R^+R^) 

= (/.(R^) + drxp(r)(r - R^).V0(R^) + 0{\/') 

= m^) + 0{V). (40) 
Performing similar Taylor expansions we obtain easily 



Also 



R-u + Rjy 



R/i + Ri- 



0(V). 



(R^) + 0(R,) 



o(v2 



and, therefore. 



0(V). 



After some algebra it is easy to show that 



Finally, 



c^. + 0(V). 



(41) 
(42) 

(43) 
(44) 

(45) 



(38) By using these Taylor approximations in Eqns. (|3q), 
(p?!), and (^8|) we obtain the final Voronoi finite volume 
discrete equations for the inviscid fluid. 



A^^ [p]^ + [p]^ 

2^ ^ c^--(M^ - [v] J, 



[^ + [^1 



■fiiy 



V- Af,^ [p]^ + [p]^ [v] + [v]^ 
^ , 2 2 '''"'■(M/^ " M-)' 



EI^^H^-.-(M.-M.)- (46) 



These equations become closed equations for M^, P^, 
by using 

r 1 ^■^A' 



(47) 



where in the last equation for the pressure we have used 
again a Taylor expansion and neglected terms of order 
0(V). 



7 



D. Finite volumes for a viscous fluid 

Having studied the inviscid fluid, in which there are 
no dissipative contributions to the motion of the fluid, 
we turn now to the general viscous fluid. The continuum 
equations are given by |^3| 

dtp{r,t)^-V-p{r,t)v{r,t), 

dtsir, t) = -VP(r, t) - V-g(r, t)v(r, i) - V- (H + HI), 
dts{r,t) = -V-s(r,t)v(r,i) 

-^V-J-' + ^W: W+|;(V.v)2, (48) 



where, by comparison with Eqns. (|T|), (32), ( |33| ) we can 
recognize the purely irreversible terms in the momentum 
and entropy equations. Here, T is the temperature field 
(which, as the pressure, is a function of p, s through the 
local equilibrium assumption). The double dot implies 
double contraction. These equations have to be supple- 
mented with the constitutive equations for the traceless 
symmetric part 11 of the viscous stress tensor, the trace 
n of the viscous stress tensor, and the heat flux J'. They 
are 



n = -27?Vv, 
H=-CV-v, 

J" 



T' 



(49) 



where the traceless symmetric part of the velocity gradi- 
ent tensor is 



Vv - i (Vv + (Vv)^) 



1 

D 



-Vv. 



(50) 



Here, D is the dimension of physical space. In principle, 
the shear viscosity rj, the bulk viscosity ( and the thermal 
conductivity k might depend on the state of the fluid 
through p, s. 

Following identical steps as for the inviscid fluid, we 
see that the viscous (irreversible) term in the momentum 
balance equation translates into 



(51) 



We now consider the gradient expansion on each term. 
For example, 

IrjVv] + \rjV^ 



Therefore, Eqn. (|l]) becomes 



(52) 



(53) 



where we have introduced 



Note that for any quantity </> we have 



(54) 



(55) 



Therefore, we see that $7^1/ is a sort of discrete version of 
the gradient operator. This discrete gradient satisfies 



(56) 



which is essentially the statement of the divergence the- 
orem as can be shown from the identity 

f d f d 

= J drx^(r) — l = dr—x^{r) 

= -Y ^M'^e^^ (57) 

where Eqn. (|2^) has been used in the last equality. 

The a, (5 component of the tensor in the first term in 
the Ihs of (BSl) becomes 



[7?Vv"\ = M,jVv"\ + 0(V) 



(58) 



In a similar way. 



[C 



(59) 



By introducing the following discrete versions of the 
quantities n, H, Vv, V-v in Eqns. (Hfl) 



n 



H„ 



[C] 



(60) 



we can write the irreversible part of the momentum equa- 
tion as 
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^f2^,-(n, + n„i) + o(v). 



(61) 



After a very similar procedure, the irreversible part of 
the dynamics in the entropy equation can be cast also in 
the form 



where 



(62) 



(63) 



Addition of the irreversible parts Eqn. (|6l|), ( |62| ) to the 
reversible part Eqns. (^6|) leads to the final Lagrangian 
finite volume discretization of continuum hydrodynam- 
ics. The numerical solution of the resulting set of ordi- 
nary differential equations would produce results which 
are accurate to order r. 

Admittedly, there are many other possibilities for ap- 
proximating equation (^l]) to first order in gradients. The 
presentation above is just a convenient form which has 
the particular GENERIC structure, as will be shown in 
the following sections. 

IV. GENERIC MODEL OF FLUID PARTICLES 

We have presented in Ref. a description of a New- 
tonian fluid in terms of discrete fluid particles. In this 
section, we present a slight modification of the fiuid par- 
ticle model presented in Ref. inspired by the results 
of the previous section on the finite Voronoi volume dis- 
cretization. 

In the fluid particle model of Ref. jl^ , the fiuid parti- 
cles are understood as thermodynamic subsystems which 
move with the flow. The state of the system was de- 
scribed by the set of variables x — {R^, P^, V^, S^, i = 
1, . . . , M}, where M is the number of fiuid particles and 
R^i is the position, is the momentum, V^i is the 
volume, and is the entropy of the fi-th fluid parti- 
cle. Because each fiuid particle is understood as a ther- 
modynamic subsystem, it has a well-defined thermody- 
namic fundamental equation. The fundamental equa- 
tion relates the internal energy £^ of the fluid particle 
with its mass M^, volume and entropy 5*^, this is 
Sfi — £iM^,Vp,, Sfj_). The local equilibrium hypothesis 
assumes that the fundamental equation for the fiuid par- 
ticles has the same functional form as the fundamental 
equation for the whole system at equilibrium. 

The volume was considered in Ref. | p5[ as an in- 
dependent variable to be included in x. In the appendix 
XI we show that, despite our intention of considering the 



volume as an independent variable, due to the particular 
form of the matrix L selected in Ref. , the volume is 
actually a function of the positions. For this reason, in 
this paper we consider from the very beginning that the 
volume Vfi of the fiuid particles is a function of the po- 
sitions of the particles, i.e., — Vp(Ri, . . . ,Ra/). We 
will actually assume that the volume of particle /x is the 
volume of the Voronoi cell of this particle, this is, Eqn. 

(lil)- 

Unfortunately, the fact that the volume is not an actual 
independent thermodynamic variable limits the applica- 
bility of the model in Ref. [|l5|. In order to recover the 
thermodynamic versatility required, it is necessary to de- 
fine a model in which the mass of the particles changes. 
From the point of view of the finite volume method of 
the previous section, it is fairly clear that the mass of 
the Voronoi cells should be taken as a dynamical vari- 
able. 

The state of the fiuid is given, therefore, by a; = 
{R^, Pp, Mfj_, Sfj,}. The energy and entropy functions are 
postulated to have the form 



S{x)^J2^^- 



(64) 



where is an implicit function of the positions of the 
fiuid particles. Regarding the dynamical invariants of the 
system, we require that the total mass M{x) = il/^ 
and total momentum P(x) = '^^P^ are the only dy- 
namical invariants. Conservation of angular momentum 
would require the introduction of spin variables in this 
discrete model 

The gradients of energy and entropy are given by 



I Z^7 OR,. ^1 \ 



dE 

dx 



\ 



dS_ 

dx 



J 









(65) 



where we have introduced the velocity, pressure, chemical 
potential per unit mass, and temperature according to 
the usual definitions, 



-P. 



p. 



(66) 
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V. REVERSIBLE DYNAMICS 

In this section we consider the reversible part ol the 
dynamics for the fluid particle model. The matrix L is 
made of M x M blocks L^^ of size 8x8. The antisym- 
metry of L translates into L^,y = —"L^^. 

We have a first strong requirements for the form of L. 
We wish that the reversible part of the dynamics pro- 
duces the following equations of motion for the positions 



R,, 



(67) 



The simplest non-trivial reversible part that produces the 
above equation has the following form 



I Z^7 OR,, ^7 \ 







p. 















where the block L^^ has the structure 



(68) 



/ 

-IS,, 



V 



16, 



fJ.1^ 



■ VfJ, 










\ 

r 



/ 



(69) 



The first row of L^j^ ensures the equation of motion (67). 
The first column is fixed by antisymmetry of L. Note 
that in order to have antisymmetry of L^j^ (which, in 
turn, ensures energy conservation), it is necessary that 
-Ai/^. Performing the matrix multiplication in 



AJ. = 
Eqn. (^ 
form 



the reversible part of the dynamics takes the 



R„ 



T T 



(70) 



We now develop the pressure term by using Eqn. (|129| ) 
of appendix |x| 



Pu 



dV, 



^ 9R 



E 

1^7^ /J 



I ^XM(r)x.(r)(r-R^)(P,,-P.) 



Pu+Pu 



+E "^^^ 



^^v{P^ Pv^ 



(71) 



where use has been made of the property 
The momentum equation becomes 



P -VA e 5^ 



P. 



E 



/ , ^Cf,4P^-P,) + Af,,fi, + r^,T,]. (72) 

Now, the basic question to answer is. What forms for 
A^y, Aui/ and F^^ should we use in order to consider 



Eqns. ([70|) as a discrete version of hydrodynamics? In 
what follows we will propose forms for these quantities in 
such a way that Eqns. ( |70| ) and (36), (|37|), ( |38| ) coincide 
as much as possible. 

The vectors A^i^ and V are easily identified by com- 
paring the mass and entropy equations in ( ^ ) and ([T^). 
The matrix A^,y is obtained by inspection from the com- 
parison between the momentum equation in ( ^ ) and 
([70|). Our proposals are 



A,, 



—c -6 



Afiv + Pv 



P^v 



Plia 2 



A^a Pfj, + Pa 
Pua 2 



(73) 



Note that — — A^^ — A^^ and, therefore, the anti- 
symmetry of L is ensured. Note also that J2y ^tJ-i' — 
and, therefore, the degeneracy condition L-dS/dx = is 
satisfied. 

By substitution of these forms into the mass and en- 
tropy equations in (|7^) one obtains the mass and en- 
tropy equations obtained in the finite volume method 
Eqns. (^6|). Substitution into the momentum equation 
( [t^ ) leads to the finite volume momentum equation (|4^), 
with an additional term which is 

^'\-Py)-P-^^{^i,~^..) 



(74) 
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This term is strongly reminiscent of the Gibbs-Duhem 
relation which, in differential forms is dP — pdfi— sdT = 0. 
For this reason, we expect that this term, although not 
exactly zero, will be very small. 

In summary, the proposed GENERIC equations for 
the reversible part of the evolution of the variables 
R/j: Pfij are 



P -VA e ^k±iii 



E 



1/ 



Rav 2 2 



J-, Cp,y ( {Pfi Pi,) 



V.) 



Afii, Pp + pi/ 



On. — 



^ -Rui/ 2 



(75) 



Here, P,j/M^, p^, = M^,/V^, and = S'^/V^. 

These GENERIC equations for the reversible part of the 
dynamics Eqns. (|7^) are identical (except for the small 
Gibbs-Duhcm term) to the finite volume discretization 
of the continuum equations of inviscid hydrodynamics 
Eqns. ( p6[ ) and can, therefore, be considered as a proper 
discretization of the continuum equations of hydrody- 
namics. Total mass, momentum, and energy are con- 
served exactly and the total entropy does not change in 
time due to this reversible motion. 



continuous description of a simple fluid can be under- 
stood with the method of projection operators. As it is 
well-known, the method produces Green-Kubo expres- 
sions for the transport coefficients. This Green-Kubo 
forms involve the correlation of the projected currents. 
Since, in the continuum case, the time derivative of the 
microscopic density field is precisely (minus) the diver- 
gence of the microscopic momentum density field (which 
is itself a relevant variable), it turns out that the pro- 
jected current vanish exactly and there is no Green-Kubo 
transport coefficient in the mass equation. The situation 
is different when one considers discrete variables. The 
discrete variables are the mass, momentum and internal 
energy of the Voronoi cells as functions of the position 
and momenta of the fluid molecules. Even though a pro- 
jection operator derivation of the equations of motion 
for these variables is extremely involved, it is possible to 
show that the projected mass current does not strictly 
vanish. This amounts to accept that the mass in a given 
cell fluctuates not only due to the indirect action of the 
random stress and heat flux but also through the direct 
effect of a random mass flux. For the time being and 
for the sake of simplicity, however, we assume that this 
random mass flux can be neglected. In this case, the 
noise term in the equation of motion (^ has the form 

(o,dP^,0,d5^). 
In the following subsections we consider two differ- 
ent implementations of the noise terms. The first one, 
through a random stress tensor and random heat flux, 
can be considered as the natural way of constructing dis- 
crete equations that, in the continuum limit, converge to- 
wards the equations of continuum hydrodynamics. The 
second implementation is a cartoon of the first one and 
leads to the Dissipative Particle Dynamics algorithm. 



VI. IRREVERSIBLE DYNAMICS 

In this section we consider the irreversible part of the 
dynamics M-dS/dx. We will postulate the random terms 
dx for the discrete equations and will construct, through 
the fluctuation-dissipation theorem (p^), the matrix M 
and the irreversible part of the dynamics. If we guess 
correctly the random terms, the resulting discrete equa- 
tions should consistently produce the correct dissipative 
part of the dynamics. 

Thermal fluctuations are introduced into the contin- 
uum equations of hydrodynamics through the divergence 
of a random stress tensor and a random heat flux p4[ , 
p5[ . In principle, one could think about a random mass 
flux that would, according to the fluctuation-dissipation 
theorem, produce an irreversible term in the mass bal- 
ance equation. Such term is absent in simple fluids but 
not in mixtures (it produces the diffusion terms). The 
reason why there is no such a random mass flux in the 



A. Finite volume hydrodynamic 

By analogy with the continuum fluctuating hydrody- 
namics we construct the random terms dP^ , dSf^ as the 
discrete divergences of a random flux 

dPp = y^^Uf_i_,y-dd-„, 

dS,^^Y.^^•^■'^~^'^^'^^^^■■T.^-^^''^■ (76) 

We will select the form (^J) for $7^^, but the particular 
form is not important for the time being. The random 
stress dcr^ and random heat flux dJ are defined by 

5 1 

do-^ = a^dW^ + 6^— tr[dW^], 
dJ« = c^dV^. (77) 
The coefficients ,b^,c^ are given by 
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1/2 



6„, = ( 2DkBT^^ 



1/2 



c„ =r„. (2fcB^ 



1/2 



(78) 



Here, D is the physical dimension of space, 77^ is the 
shear viscosity, is the bulk viscosity, and is the 
thermal conductivity. These transport coefficient might 
depend in general on the thermodynamic state of the 
fluid particle /z. The particular form of the coefficients in 
Eqn. ( [tsI) might appear somehow arbitrary. Actually, it 
is only after writing up the final discrete equations and 
comparing them with the finite volume equations (61), 
( |62|) that we could extract the particular functional form 
of these coefficients. 

The traceless symmetric random matrix dW ^ is given 

by 

^l^l ['^W^ + '^W^] - ;^tr[rfW^]l. (79) 



dW^ is a matrix of independent Wiener increments. The 
vector dVf^ is also a vector of independent Wiener incre- 
ments. They satisfy the Ito mnemotechnical rules 



0, 



(80) 



where latin indices denote tensorial components. Note 
that the postulated forms for dP^ , dS^ in Eqn. ( |76| ) sat- 
isfy 



Y^Vf,-dPf, + T^dS,, = 0, 

Y.dP^ = o, 



(81) 



and, therefore, Eqns. ( pi] ) are satisfied. This means that 
the postulated noise terms conserve momentum and en- 
ergy exactly. It is now a matter of algebra to construct 
the dyadic dxdSp^ and from Eqn. ( [l0| ) extract the matrix 
M. The procedure is rather cumbersome but standard. 

Once M is constructed, the terms M-dS/dx in the 
equation of motion (|l|) can be written up. By assum- 
ing that the transport coefficients do not depend on the 
entropy density (but they might depend on the mass den- 
sity), the resulting equations of motion are 



Vi 



dt 



kBTf_, 



Tf^dSf,. 



2D 



(82) 



In these equations, we have introduced the same quanti- 
ties as in Eqn. (|60|). The heat flux J-^ is defined by 



(83) 



Finally, the heat capacity at constant volume of particle 
/i is defined by 



T,. I S- 



(84) 



We observe that, quite remarkably, the above equa- 
tions are in the limit fc^ ^ identical to the irreversible 
part of the particular finite volume discretization of the 
continuum hydrodynamic equations presented in section 



III 



We have, therefore, shown that these equations ( |82| ) 
are a proper discretization of the irreversible part of hy- 
drodynamics with thermal noise included consistently. 

By collecting the reversible part L-dE/dx in Eqn. ( |6g 
and the irreversible part M-dS/dx in Eqn. (^2|), the 
final equations of motion for the discrete hydrodynamic 
variables could be finally written. 



B. Irreversible part of DPD 

In this section we show how, by postulating a different 
form for the noise dx^ one can obtain an irreversible part 
which is closely related to the irreversible part of the 
Dissipative Particle Dynamics model. The Dissipative 
Particle Dynamics model that we present here, then, is a 
natural generalization of the classical DPD model 0] in 
which not only an internal energy (or entropy) variable 
is included as in Refs. 1^ but also a mass density variable 
is introduced. From the GENERIC point of view, this 
DPD mod el diff ers from the finite volume hydrodynamics 
in section VI A model only in the form of the dissipative 
and random terms. 

Instead of ( [76| ) the postulated structure of the random 
terms is dx (0, dP^, 0, dS"^) with the following defini- 
tions 



dP^ =^B^,dW^„ 



dS^ = -2^5I^A'--^A'-^^M- + T^^^M-rfl^.., (85) 
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where B,, 



0, and Aai, 



A,. 



suitable functions of position and, perhaps, other state 
variables. The independent Wiener processes satisfy 
dW^i, = dW^^ and dV^^, = ~dVuf_i and the following Ito 
mnemotechnical rules 



dWf_,^>dWuv' = [5^^u5^i'u' + Sf^„>df_,>t,]dt, 
dW^^,dV,,' =0. 



(86) 



Note that the noise terms in Eqn. ( p5| ) satisfy the re- 
quirements (|ll|) which take the form 



(87) 



The first equation ensures momentum conservation while 
the second equation ensures energy conservation. Note 
that the random force dP^ provides "kicks" to the par- 
ticles along the line joining the particles, and it satis- 
fies Newton's third law. The first term of the random 
term dSf^ is suggested by the last equation in Eqn. ( |87| ) 
whereas the last term is dictated by our wish of modeling 
heat conduction Q. 

We note that the DPD noise terms ( ^5[ ) can be viewed 
as a cartoon of the noise terms ( [76[ ) in the finite Voronoi 
volume model. For this reason, we will assume the fol- 
lowing form for A^^, B^j^ to remain as close as possible to 
{ ffdj } while retaining the correct symmetries for A^i,, 'B^^, 



V 



1/2 



B, 



1/2 



2kB ^''^^ 1 1 n^,. 



(88) 



We have introduced parameter with dimensions of 

a thermal conductivity and a coefficient 7 with dimen- 
sions of a viscosity. This DPD model has only a single 
viscosity, instead of two viscosities (shear and bulk) ap- 
pearing in the finite volume model of the previous sub- 
section. We denote with V = Vr/M the average volume 
per particle. 

It is now a matter of simple algebra to construct 
the irreversible matrix M for the DPD algorithm with 
the fluctuation-dissipation theorem (|l0|). The final irre- 
versible part of the equations of motion are 

T,dS,\,^.^ =^^(0^..v,.)2di 

-i-^^f|2^(r.-r^)dt 



mV 



J2^l.T,dt + Tf,dSf,. 



we have defined the pair viscosity as 



7^1^ = 7 1 - 



ks 



(89) 



(90) 



We discuss each term of these equations now. The mo- 
mentum of a particle changes irreversibly due to a fric- 
tion force that depends on the velocity differences be- 
tween particles and due to a random noise. This is the 
conventional form of the DPD forces except for the fact 
that the friction coefficient depends on the temperatures 
of the particles (although with a small prefactor of order 
fc-B /C'ai) ■ Concerning the equation for the evolution of the 
entropy, the first term models the process of viscous heat- 
ing, the fact that the motion of the particles creates an 
internal friction that increases the internal energy of the 
particles. This term is dictated basically by the require- 
ment of energy conservation. The second term takes into 
account the process of heat conduction and it can be un- 
derstood as a simple discretization of the heat conduction 
equation. Note again that this form of the conduction 
term is physically more reasonable that the expressions 
given in |^]. The next three terms are proportional to fc^ 
and ensure the exact conservation of energy, as can be 
seen by explicitly computing dE = ^dx + ^-^^^dxdx. 

At this point, we would like to make a close compar- 
ison between the model presented in Ref. and the 
one presented in this paper. In Ref. |Q the variables 
used to describe the system of fluid particles are the po- 
sitions, momentum and energy of the particles. The mass 
is assumed to be a constant. We will focuse on the mo- 
mentum equation because in order to compare the energy 
equation in both works we should transform our system 
of variables with the entropy to a system with the en- 
ergy as independent variable. The momentum equation 
in 1 23 1 is (changing their notation to our notation) 



2 



^ [vp^ + (v^^ •e^^)e^^] \ + M^g + (91) 



where g is a body force like gravity and is a stochastic 
force with a form fixed by the fluctuation-dissipation the- 
orem. We note that the reversible part of this equation, 
the pressure term, is identical to ours except for the small 
Gibbs-Duhem term required for thermodynamic consis- 
tency when the mass and entropy are not conserved by 
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the reversible part of the dynamics as happens in our 
case. The pressure difference instead of the sum is used 
in Eqn. (|9^), but both are equivalent in view of Eqn. 
(p7|). The irreversible term proportional to the "viscos- 
ity" 77, is of a form similar to the DPD irreversible term 



R„ 



dP^lii-r in Eqns. (g9|), with fi^^, given by Eqn. (|5j). A 
term proportional to v^i, , which breaks angular momen- 
tum conservation is included |^ but poses no conceptual 
problems. However, we note that dimensionally rj is not 
a true viscosity, a volume factor is missing. Actually, we 
expect that a simulation of Eqns. ( |9l| ) with a fixed value 
of 1] at different resolutions (this is, the same external 
length scale, channel width for example, and different 
number of fluid particles which have, consequently, dif- 
ferent typical volumes) would produce different values for 
the actual viscosity of the fluid being modeled. 



VII. THE SPH MODEL 

In this section we show that the Smoothed Particle 
Hydrodynamics algorithm for an inviscid fluid is a par- 
ticular case of the GENERIC Eqs. ([70|). For a viscous 
fluid, the SPH algorithm involves the irreversible part 
( ^ ) with no fluctuation effects, i.e. fc^ = 0. 

In SPH, a density variable associated to each fluid 
particle is defined through 



(92) 



where the weight function W{r) is a bell-shaped function 
with finite support h and normalized to unity. If particles 
are close together, then the density defined by Eqn. (^) 
is higher in that region of space. From the density one 
can define a volume associated to each fluid particle 



Vp(Ri, . . . ,Rm) — -5- 



(93) 



The SPH algorithm for an inviscid fluid is obtained from 
Eqns. ( [70I) by using A^^= A^i, — T^i, = 0. The deriva- 
tive of the volume (|9^) with respect to the positions of 
the particles is easily computed as 



dl 



(94) 



M. = 0, 



(96) 



This is the form of the discretization preferred by Mon- 
agan jlj] . It is remarkable that this form is forced solely 
by the selection of the volume function in Eqn. ( |93| ) and 
the GENERIC structure. For a viscous fluid, the irre- 
versible part of the SPH equations are simply obtained 
from Eqns. ( ^4) b y using as fJ^y the function —dVv/d'R.^ 
given in Eqn. (|94|). 

In the SPH equations (|96| ) , the mass o f the particles is 
a constant. We will show in section VIII, when we discuss 



the equilibrium distribution function, that this constancy 
makes the algorithm unsuitable for studying gas-liquid 
coexistence. On the other hand, the SPH algorithm suf- 
fers from an unphysical feature: if the pressures of the 
neighbors of particle are equal to the pressure of this 
particle /i, there is still a remnant force on this particle, 
which is physically unacceptable. A possible correction 
of the above defect is as follows. 



Instead of the volume defined as in (93) we define the 
volume as 



Vt 



(97) 



which satisfies '^^V^ = Vt- Note that the correction 
factor Vt/ X]fc ^fc ^ is expected to be close to 1. Note also 
that now the volume of a given particle depends on the 
coordinates of all the particles in the system. This breaks 
the local definition of the volume of a particle in terms 
of the positions of the neighbours. 

Now, by using this expression for the corrected volume 
of a particle we obtain 



where 



and, therefore, 



iP,-P)^ + {P.-P)^ 

dn dy 



(98) 



a;^., (99) 



VF'(r^^)e^j,. 



(95) 



The prime here denotes the derivative and e^^ is the unit 
vector joining particles 

Substitution of Eqn. ( p^ ) into the reversible equations 
(^) leads to 



where P is a sort of "spatial average" of the pressures of 
the particles, i.e.. 



(100) 



Equation (£^) has the following good features: Total mo- 
mentum and total volume are conserved variables. If the 
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pressure of all the particles is exactly the same, there 
is no force on the particles. And the following draw- 
backs: The force on particle fji depends on the state of 
all the particles of the systems. Therefore, the change in 
momentum of a bulk of fluid particles takes place non- 
locally, not through the neighbourhood of this bulk. This 
non-local effect breaks the local transport of momentum 
and therefore the macroscopic hydrodynamic behaviour. 
We expect that this effect is small, particularly when the 
range of the weight function is much larger than the in- 
terparticle distance. The second drawback is that the 
remnant force is zero only if the pressures of all the par- 
ticles of the system are equal. If in a local region of the 
system the pressures are equal for the particles in that 
region but different to the pressure of other regions, then 
there still exists a remnant force on the particles of that 
region. 

A final word on the overlapping coefficient is in or- 
der. We define the overlapping cocfhcient s as the ratio 
of the range h of the weight function W{r) to the typ- 
ical interparticle distance A between fluid particles, this 
is s = h/X. When s 3> 1 a given particle has typically 
many neigbhours. It is clear that if the overlapping coef- 
ficient is much larger than 1 then, for homogeneously dis- 
tributed particles all the volumes of all the particles will 
be typically the same. In the limit s large, then, the pres- 
sures of the particles, which depend only on the volumes 
of the particles, will be the same and the forces on the 
particles will be negligible. In this limit, the equations of 
motion ( p6| ) lose their sense. Actually, we expect that the 
volume defined through ( p7| ) will have sense only if the 
overlapping coefhcient is slightly larger than 1. However, 
the usual derivations of SPH by convoluting the contin- 
uum equations of hydrodynamics with the weight func- 
tion W{r) make sense only in the limit of large overlap- 
ping, in such a way that the integrals can be reasonably 
approximated by sums. Because of this inconsistency and 
the problems discussed above, wc prefer the formulation 
of discrete hydrodynamics in terms of finite Voronoi vol- 
umes rather than in terms of weight functions as it is 
done in SPH. From a computational point of view, we 
note that the large cpu time required in order to update 
the Voronoi mesh (i.e. compute the quantities A^j^^^c^^) 
may be compensated by the fact that typically only six 
neighbours need to be considered in the finite Voronoi 
volume simulation whereas 30-40 neighbours are needed 
in a SPH simulation in 2D. 

We would like to comment finally on the approach pro- 
posed in Ref. where a finite volume approach similar 
to the one presented here is advocated. However, these 
authors do not consider the Voronoi limit a Q but 
rather take a finite width for the function A(r), provid- 
ing an overlapping coefhcient of typically 1.4. Because 
they do not consider the Voronoi construction, they en- 
counter the difficulty of evaluating integrals similar to 
those in Eqns. In one dimension, as is the case con- 



sidered in Ref. , this can be achieved with a numerical 
integration method, but this becomes readily infeasible 
in higher dimensions. 



VIII. EQUILIBRIUM DISTRIBUTION 
FUNCTION 



In this section we discuss the equilibrium distribution 
function p°'^(x) corresponding to the equations of mo- 
tion ( [75[ ) plus the irreversible terms (p2[). Note that the 
equilibrium distribution function of Eqns. ( ^2|) and (|89| ) 
is the same, irrespective of the actual form of the irre- 
versible part of the dynamics. This is because both sets of 
equations have the GENERIC structure. The GENERIC 
structure of the equations of motion ensures that the 
equilibrium distribution function for these variables is 
given by Einstein distribution function in the presence 
of dynamical invariants, Eqn. (^. Because total mass 
M {x) total energy Eq and momentum Pq are conserved 
by the dynamics, the Einstein distribution function will 
be given by 

p^\x) = ^S{M{x) - Mo)6{E{x) - Eo)6{P{x) - Pq), 



X e-xp{kg^S{x)}, 



(101) 



where we have assumed that we know with absolute pre- 
cision the values of the total mass Ado, energy Eq and 
momentum Pq at the initial time. This is the situation 
in a computer simulation, is a normalization factor 
that ensures the normalization of p''^{x). A word is in 
order about the total volume. We note that due to the 
definition of the volume in Eqn. (|l5|), any configuration 
of positions of the particles gives that the total vol- 
ume has the same value ~ ^o- Therefore, even 
though total volume is conserved, it does not produce a 
restriction in the form of a delta function in Eqn. (101). 



A. The most probable state 

The most probable state x* at equilibrium according 
to ( |l0l[ ) is the one that maximizes the entropy S{x) sub- 
ject to the constraints AI{x) = Mq, E{x) = Eq, and 
P(a;) = Pq. By introducing Lagrange multipliers /3, A 
and V, the most probable state x* is the state that max- 
imizes kg^S{x) - (3{E{x) - V-T'{x) - XM{x)) without 
constraints. By equating the partial derivatives with 
respect to every variable to zero, one obtains the fol- 
lowing implicit equations for the most probable values 

i' = {b-;,p;.m;,s;) 
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-I- = v 

M* 



(102) 



The second equation states that in the most probable 
state all particles move at the same velocity V which 
might be set to zero without loss of generality. The 
two last equations state, then, that the temperature and 
chemical potential per unit mass of all the fluid parti- 
cles are equal at the most probable value of the discrete 
hydrodynamic variables. This implies that the pressure 
is also the same for all the fluid particles (in a simple 
fluid the intensive parameters are not independent [p6|). 
The first equation is, therefore, trivially satisfied (be- 
cause Vu = ctn). 



B. Marginal distribution functions 

In this subsection we will integrate out the momen- 
tum variables in Eqn. (101) in order to have more spe- 
cific information about the distribution of different vari- 
ables at equilibrium. To this end, we denote the state 
X = (y, {P}) where y = ({R}, {Af}, {S}) is the set of po- 
sitions, masses, and entropies of all particles. Note that 
the total entropy and internal energy in Eqn. ( p^ ) do not 
depend on momentum variables. 

By integrating the distribution function p'^'i{x) over 
momenta we will have the probability p'"^{y) of a real- 
ization of y 



p^\y) = exp {S{y)/kB} ^5 - 



X / dPi 



M 



MS 



2AL 



+ £^{y) -Eo 



exp {S{y)/kB} ^5 - Mo^ 



M 



2 



£;o-^f^(2/) 



where we have used Eqn. (153) of appendix KIl 



We find now a convenient approximation to Eqn. (103) 
by noting that this probability is expected to be highly 
peaked around the most probable state. Therefore, for 
those values of the variables £fi{y) for which p'^^(y) is 
appreciably different from zero we can approximate 



-I p 



-l p 



exp{/3*5](£;-£^)} 



ctn. exp{-/3*y^g^}. 



(104) 



where £* is the most probable value of and P = 
D{M — l)/2 — lisa very large number and we have 
introduced 



(3* = 



i:i(Af-i)/2-i 



DM/2 



Finally, we can write Eqn. ( |l03|) as 



(105) 



iy)^^exp^^S{y)/kB~P*J2^^^(y^ 



(106) 



where fl'f^ i s the corresponding normalization function. In 
Eqn. ( |l06|) we have neglected a term log A/^ in front 
of J2fi ^/j(y)- Note that is a first order function of its 
arguments S"^, Af^, and, therefore, it is of order A/^. 

Changing back to the original notation we write Eqn. 
(106) in the form 



p<^'i({R, Af, S}) = exp 1^ S^/kB ~ P*£^{M^, S^, V^)| 

X 5 - . (107; 

We see, therefore, that by integrating the momenta the 
"microcanonical" form Eqn. (101) becomes the "canoni- 
cal" form ( p^ ). 

We are interested now on the distribution function 
P{M^, Sfi) that the particular cell p, has the values 
Af^, for its mass and entropy, irrespective of the values 
of th e rest of the variables in the system. We integrate 
( 107 ) over the variables of all cells except Af^ , 5^ . 



(103) 



X exp 1^ - /3*f.(Af„ S,, V.)| 

X 5 (j2m^-mA 
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exp 1^ 5./fcs - f3*SAMu, K)| 



M 



■■■,Vm) 



where we have introduced the identity 

„ M 
J d{V}Y[S{V,-V,i{K}))^l, 



(108) 



(109) 



and the function 

1 r 

F{Vi,...,Vm)^ttm / d{Jl}l[S{V,-V,{{R})). 

(110) 

This function is the probabihty density that the the par- 
ticles have the particular distribution V^i , . . . , Vm of vol- 
umes provided that the distribution function of the po- 
sitions is uniform. The calculation of this function is 
difficult and we do not attempt to do it. Rather, we 
will assume that this function is highly peaked around 
{V , ■■■,¥) in such a way that all the cells have approxi- 
mately the same volume V = Vt/M. Under this approx- 
imation, Eqn. (108) becomes 



P(M^, S^) = exp {S^^/kB - /3*fp(M^, 5^, V)} 



X $(Mo - M^), 
where we have introduced the function 

M 



(111) 



exp <^ J2 ^'^I^B ~ P*£u{M,,S,, V) 




(112) 

The functional form of ^{X) is very well approximated 
by an exponential. This can be seen by taking the deriva- 
tive 



$'(X) = - ^ d(*^-i){M}d(^-i){^} 



M 



exp <j SJkB - I3*£AM,, S,, V) 



X) 



= -P* / d^^^-^^{M}d^^^-^^{S}iiiAU,S^,V) 



M 



exp <j Y ^-/^s - P*SAM,, S,, V) 



M 

xS\J2Mu- X) 



(113) 

where we have integrated by parts in the second equality 
and a is the label of any cell exc ept fi . The most probable 
value of the integrand in Eqn. ( |ll3| ) is the solution of 



(114) 



where the chemical potential per unit mass equates A, a 
suitable Lagrange multiplier that accounts for the mass 
conserving delta function. We expect that when the num- 
ber of variables M is very large, the integrand in Eqn. 
(113) becomes highly peaked around this most probable 
value. In this case, we have 



$'(X) « -(3*\^{X), 



(115) 



and therefore ^{X) is the e xpo nential exp{— /3*AX}. 
Returning back to Eqn. (HI) we have finally. 



P(M^, S^) = - exp{S^/kB - 13* {£{M^, V) - XAQ } , 

(116) 

where Z is the appropriate normalization. Note that the 
argument of the exponential is the Gibbs free energy with 
fixed values for the inverse temperature (3* and chemical 
potential per unit mass A. 



C. van der Waals fluid 

We will particularize the discussion of P{M^,S^) by 
considering the fundamental equation £{Mf^, Sf^^V) for 
a van der Waals fluid. First we note that the internal 
energy is a first order function of its variables and, there- 
fore. 



£{M^,S^,V)^Ve{n^,s^), 



(117) 



where — M^/{niQV) is the number density, — 
Sfi/V is the entropy density, and e = £/V is the inter- 



nal energy density of cell /i. Therefore, from Eqn. (llC) 
we can obtain the probability density that cell fi has the 
values n, s for its number density and entropy density. It 
is given by 
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P(n,s) = lexpy {s/fcs -/3*(e(n,s) - Xmon)} . (118) 
Zj 

It is convenient to use reduced units for the van der 



Waals fluid (see appendix XIII for details of notation). 
In reduced units the distribution function becomes 

P(n, S) = ^ exp{F \l - /3'=''*(?(n, ~s) - fi^^'^n 

The most probable value of this distribution function 
occurs at n* , s* which are solutions of the equations 



)} (119) 



dfi 
de{h*,s*) 
ds 



jl{n* , s* 



~cxt 



1 



= T{n*~s*) = — 



(120) 



For fixed ri, the relation between T and s is monotonic, 
whereas the chemical potential has the form 



FIG. 4. Zoom of the previous figure for the isotherm 
T = 0.85. The equal area construction gives a value 
^cxt _ _i4^ ggYl which produces two humps of equal height 
in P{n,f). 



In Fig. H we show the chemical potential for different 
values of the temperature. We observe that for T < 1 
there are three solutions for n* in the first equation in 
(120). This means that depending on and the 



distribution function s) can present a bimodal form. 

Due to the particular functional form of the fundamen- 
tal equation for the van der Waals gas, it is more conve- 
nient to study the distribution function P{n, T) instead 
of P(h,s). This function is computed in the appendix 

xnt 



/i = T In 



3-fiJ 4 2 I c 



(121) 




FIG. 3. Chemical potential fl{n, T) as a function of n for 
difi'erent T. In descending order, T = 0.85, 1.0, 1.15. Observe 
that for T < 1 the equation /i(n, T) = /i™* might have three 
solutions for n, depending on the actual value of /i™*. 



P (n,T) 




2.5 



FIG. 5. The distribution function P{h, T) for an external 
temperature = 0.85. The external chemical potential 

is /i^"' = -14.85. The typical volume of the cells is = 20. 
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P (n,T) 




2 . 5 



FIG. 6. The distribution function P{n, T) for an external 
temperature 1/(3'^^^ = 0.85. The external chemical potential 
is /i""* = -14.8971. The typical volume of the cells isV ^ 20. 
The two maxima have equal height. 




2 . 5 



sition, because a cell of size V can have a non-vanishing 
probability of having two different values of the density. 
The external chemical potential jl'^^'^ controls the relative 
magnitude of the two maxima. For the particular value 
^oxt _ —14.8971, at this value of the external tempera- 
ture, the two heights are equal. This value is called the 
chemical potential of coexistence and satisfies the equal 
area rule, see Fig. ^. 

The size V of the typical volume of a cell controls the 
sharpness of the distribution function. In Fig. ||, ^, |l^ we 
show the distribution function P(n, T = 0.85) for differ- 
ent values of V . The fluctuations become much smaller 
as V becomes larger, consistent with the fact that in the 
thermodynamic limit fluctuations vanish. Another inter- 
esting property of the thermodynamic limit is that the 
relative height of the to peaks in the distribution func- 
tion of the density increases as ^ oo. Eventually only 
one of the two peaks survives. This is true, whenever the 
two peaks are not equal in height. When they are exactly 
equal both spikes coexist in the thermodynamic limit. 



P(n,f) 




FIG. 8. The distribution function P{n, T) along the 
isotherm T = 0.85 for two different values oi V = 2Q (solid 
line) and 60 (dotted line) . Here, the external chemical poten- 
tial is /i™* = —14.85. Observe that the distribution becomes 
more peaked as the typical volume V of the cells increases. 
The relative height increases and, eventually, in the thermo- 
dynamic limit V — > 00 only one peak survives. 



FIG. 7. The distribution function P{n, T) for an external 
temperature 1/(3'^^^ = 0.85. The external chemical potential 
is /I^"" = -14.95. The typical volume of the cells is l7 = 20. 



In Figs. ^,^,0 we show the distribution of number den- 
sity and temperature P(n, T) for a value of the external 
temperature l//3°''* = 0.85 (below the unit critical tem- 
perature) and three different values of the external chem- 
ical potential jT''^ = -14.85,-14.8971,-14.95. We ob- 
serve the presence of two humps at the same value of the 
temperature T — 0.85 and different values of the density 
n. The existence of a bimodal structure in the distribu- 
tion function P(n, s) is a reflection of the gas-liquid tran- 
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0.5 1 1.5 2 2.5 3 

h 

FIG. 9. The distribution function P{n, T) along the 
isotherm f = 0.8971 for two different values of F = 20 (solid 
line) and 60 (dotted line). The chemical potential is that of 
coexistence, /i™* = —14.8971. Because the heights are equal, 
as the typical volume V of the cells increases, they remain 
equal at later times. In the thermodynamic limit two spikes 
localized at the gas and liquid densities remain. 
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FIG. 10. The distribution function P(n,f) at the crit- 
ical point. The critical isotherm T = 1 is plotted for 
^oxt ^ _i7 35 foj. two different values oi V = 2Q (solid line) 
and 60 (dotted line). Note the broad distribution of densities, 
which remain broad even in the thermodynamic limit. 



IX. DISCUSSION 

We have presented a Lagrangian finite Voronoi vol- 
ume discretization of continuum hydrodynamics and 
have shown that the discrete equations have almost the 
GENERIC structure. The GENERIC structure can be 
restored by simply adding a small term into the reversible 
part of the dynamics. Therefore, the obtained equations 
conserve mass, momentum, and energy, and the entropy 
is an strictly increasing function of time in the absence 
of fluctuations. Thermal fluctuations are consistently in- 
cluded which lead to the strict increase of the entropy 



functional and to the correct Einstein distribution func- 
tion. The size of the fluctuations is given by the typical 
size of the volumes of the particles, arguably scaling like 
the square root of this volume. The need of incorporat- 
ing thermal fluctuations in a particular system will be 
determined by the external length scales that need to be 
resolved. For example, if sub-micron colloidal particles 
are considered, we need to resolve the size of the col- 
loidal particle with fluid particles of size, say, an order of 
magnitude or two smaller than the diameter of the col- 
loidal particle. For these small volumes, fluctuations are 
important and lead to the Brownian motion of the par- 
ticle. A ping-pong ball, on the other hand requires fluid 
particles much larger, for which thermal fluctuations are 
negligible. Of course, one could use a very large number 
of small fluid particles to deal with the ping-pong ball, 
but in this case the (large) thermal fluctuations on each 
fluid particle average out among the (large) number of 
fluid particles. 

The original formulations of Dissipative Particle Dy- 
namics lack this effect of switching-off thermal fluctua- 
tions depending on the size of the fluid particles. This 
is due to the fact that early formulations did not include 
the volume and/or the mass of the particles as a relevant 
dynamical variable. In this paper, we have extended the 
range of variables and have shown that a thermodynami- 
cally consistent Dissipative Particle Dynamics model can 
be formulated. This model has the same reversible part 
as the finite Voronoi volume model. Hence, the usual 
conservative forces between dissipative particles are sub- 
stituted by truly pressure forces. Therefore, prescribed 
thermodynamics can be given, without the need of "re- 
verse inference" (find out which conservative force will 
produce the desired equation of state Q). Even though 
the present DPD model has clear advantages over previ- 
ous DPD models, it is yet inferior to the finite Voronoi 
volume model presented in this paper. Note that the fi- 
nite volume is closely related to a discretization of Navier- 
Stokes equations, whereas the DPD model in this paper 
does not converge to the Navier-Stokes equations, even 
though it displays correct hydrodynamic behaviour. This 
means that in order to simulate a fluid of a given viscos- 
ity, one has to tune the viscosity parameter of the DPD 
model by previous simulation runs. With the finite vol- 
ume model the true viscosity is given as input. 

We have shown also that the Smoothed Particle Dy- 
namics method has the GENERIC structure. However, 
we have pointed out several problems that favour the use 
of the finite volume method of the present paper. For 
example, in the Smoothed Particle Dynamics method we 
have noted that remnant forces appear even for the sit- 
uation in which all particles are at rest with identical 
pressures. This effect can be made arbitrarily small by 
increasing the overlapping coefficient. But this has two 
drawbacks. The first one is of practical nature: when the 
overlapping coefficient is increased the number of inter- 
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acting neighbours increases and, therefore, the compu- 
tational time increases too. A second, more conceptual 
drawback is that in the large overlapping limit which, 
strictly speaking, is needed in order to "deduce" the 
SPH equations for Navier-Stokes, the volume (or den- 
sity) of the particles becomes a constant and the equa- 
tions start loosing its sense. The pressures, which de- 
pend essentially on the volumes, become a constant and 
the resulting equations cannot be claimed to be a faithful 
discretization of Navier-Stokes equations. Even though 
hydrodynamic behaviour is displayed (momentum is con- 
served locally) we expect that the thermodynamics and 
the transport properties of the simulated fluid differ from 
the actual desired values. 

Concerning the GENERIC finite volume model of fluid 
particles presented in this paper, one of the most sur- 
prising realizations has been the fact that the mass of 
the fluid particles changes due to the reversible part of 
the dynamic. In Refs. ||l^, this was overlooked. Of 
course, we could, in principle, formulate a model with 



c^u = in the final discrete equations ([75|). This would 
result in a great simplification of the equations without 
losing the desired GENERIC properties. This approxi- 
mation would imply that the mass of the particles are 
constant. This spoils the equilibrium distribution func- 
tion that would not be given by the Einstein distribution 
function (101) but rather by 



M 

(^) = ^U^iM.ix) - Ml)5{E{x) - E,)5{V{x) - Po), 



exp{fc^^S'(a;)} 



(122) 



where is the initial value of the mass of particle ^. 
In this case, not only total mass but also the individ- 
ual masses are dynamical invariants of the discrete equa- 
tions when c„i, = 0. We note that we could not follow 



the same arguments that lead to ( [116D and which de- 
scribe the liquid-gas coexistence of a van der Waals fluid. 
This very same argument can be applied to the case of 
the SPH model discussed in section VIl which renders 



the SPH model unsuitable to discuss gas-liquid coexis- 
tence. From a theoretical point of view, the possibility of 
simulating liquid-gas coexistence with a particle model 
represents a stringent condition on the thermodynamic 
consistency of the model proposed. Actually, we have 
learnt that the mass, and not the volume, of the fluid 
particles should be included as an independent dynamic 
variable if coexistence is to be described. 

Finally, we would like to note that the model presente 
does not exibit the phenomena of surface tension. Surface 
tension can be easily included as an extra conservative 
contribution to the energy function [5^ . 
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X. APPENDIX: VORONOI PROPERTIES 

In this appendix we explicitly compute the derivative 
of the volume of cell fi with respect to the position Rj^ 
of cell ly, this is 

G^. = = 4 / dRx^{R){S^, - x.(R))(R - R.). 

City (7 Jy^ 

(123) 

Its worth considering the cases ^ v and fi = v explic- 
itly. 

G^, = -\( dRx;.(R)x.(R)(R-R.), v^^M 

\( dRx^(R)(i-XM(R))(R-RM) 
I rfRx^(R)x.(R)(R - R^) 



G 



(124) 



It is convenient to rewrite Eqn. (124) for i> =/= fi as 



G^. = dRx^(R)x4R) (^R- 52i±i^ 

dRx4R)XM(R) (125) 



>Vt 
1 



2f72 



with Rp,y = R^ — Rjy. The first task is to compute the 
limit cr ^ for these two integrals. For this reason, 
it is instructive to work out the actual forms of XmI-^) 
and X/j(R-)Xi'(R-) for the case that only two particles are 
present in the system, as has been done by Flekkoy and 
Coveney [|3| . Simple algebra leads to 



Xm(R)x.(R) = 



1 



1 + exp{-R^^-(R-- (R^ 
1 



R.)/2)/a2} 



4cosh^(R^^-(R- (R^ + R^)/2)/2r72}) 



(126) 



Note that X/i(R-)Xi'(R) is different from zero only around 
the boundary of the Voronoi cells of particles /i, v. In the 
limit of small a this is even more true. The integrals in 
(125) therefore can be performed not over the full volume 
Vt but only over a region df^^, "around" the boundary of 
the /i, V cells. In this region, we can further substitute the 
expression of X/i(r)x!^(R')! which depend on the positions 
of all the center cells, by Eqn. (126), which depends only 
on the position of the centers of cells /i, v. Actually, we 
can make a translation from R to R' = R — (R^-|-R^)/2 
(we put the origin exactly at the boundary between cells). 
We can also make a rotation in such a way that the x axis 



22 



is along the line joining the cell centers. In this way, we 
can write 



^ / dRx^(R)x^'(R) 

O" JVt 



1 

4^ 



dRl- 



1 



—A r 



cosh^(R'-R^^/2cr2) 
1 



dx- 



cosh^(xi?^^/2(T2) 



(127) 



Note that Jjj°° cosh '^{x)dx = 1. Here, A^y is the actual 
area of the boundary between Voronoi cells of particles 

In a similar way, one computes the first integral in Eqn. 
(125) with the result 



\ [ dRxM(R)xA.(R) (r 

^ JVt V 



R,, + R„ 



Aau 

— —i 



(128) 



where the vector c^^ is, by definition, the position of 
the center of mass of the face between Voronoi cells fi, v 
with respect to the point (R^ + R^)/2. Collecting Eqns. 
(p^),(p^ and (p^ one finally obtains 



Gnu 



(129) 



where 

























(130) 



due to Eqn. (124). Also 



Note that E^G^i' 

~ because of Eqn. ( |l23| ) and the fact that 
total volume is a constant, independent of positions. 



XI. APPENDIX: DEPENDENT VS. 
INDEPENDENT VARIABLES IN GENERIC 

What happens if in the description of the state of a 
system one introduces, perhaps without knowing it, vari- 
ables which are not independent of each other? We pro- 
vide in this appendix the answer to this question. 

Let us assume that the state of a system is described 
by a set of independent variables x. Assume, for the 
sake of simplicity, that the dynamics is purely reversible. 
The case M ^ follows along a similar way as the one 
presented here for the reversible part. 

Now assume that the basic building blocks of 
GENERIC have the following structure 



E{x)=E{x,y{x)) 
S{x) = S{x,y{x)) 
L{x) = L{x,y{x)) 



(131) 



where y{x) is a prescribed (vector) function of the sate 
X. From the chain rule one has 

VE{x) = V,E(x, y{x)) + J^{x)\/yE{x, y{x)) (132) 

and similarly for S{x) and L{x). Here, the matrix J{x) 
is defined by J(x) — dy/dx. 

The GENERIC reversible part of the dynamics is given 

by 

x = L{x)VE{x) 
= L{x,y{x))\/j,E{x,y{x)) 
+ L{x,y{x))J^{x)VyE{x,y{x)) (133) 

We can consider the time derivative of the functions 
y{x(t)) which, again through the chain rule, is given by 



y{x) = J{x)x 



(134) 



By using now Eqn. ( |133| ) into ( |134| ) we can finally group 
both equation into the form 



L LJ' 



v,£; 



y 



JL JLJ'^ I \ VyE 



(135) 



or, by introducing the obvious notation z = {x, y{x)), 

z = C{z)V,E{z) (136) 

A very striking nicety is that the matrix C{z) has all 
the required properties for being a proper GENERIC re- 
versible matrix. This is 

C^iz) = -dz) 
C{z)\7^S{z) = 

v,[/:(z)v,s(z)] = 



(137) 



as can be easily checked from the properties of L{x). 

Now, let us assume that we start describing the state 
of a system by a vector z — (x, y) and that the matrix 
L{z) has the structure 



L{x,y) L{x,y)j'^{x) 
J{x)L(x,y) J{x)L{x,y)j'^ (x) 



(138) 



with J{x) — dh/dx for some set of functions h{x). Then, 
we will have 



d 9h d 



(139) 



and, therefore, y = h{x). Therefore, if in a given model 
the matrix L(z) happens to have the structure given in 



23 



(138), then, necessarily, the variables y are dependent on 
the variables x. 

In our formulation in Ref. [T^ we have assumed that 
the volume evolves according to 



M 



(140) 



which for most reasonable forms of Vtui, has the form of 



Eqn. (139). Therefore, the dynamical equation for the 
volume is nothing else than application of the chain rule 
and the volume cannot be considered as a truly indepen- 
dent variable of positions. 



XII. APPENDIX: MOLECULAR ENSEMBLE 

In this appendix we want to compute explicitly the 
following integral 



So 




M / M \ 



(141) 



which appears repeatedly when computing molecular av- 
erages. The equation (2M^)^/^P^ = Pq are actually 
D equations (one for each component of the momentum) 
which define D planes in R^^^ . The integral in (141) 
is actually over a submanifold which is the intersection 
of the D planes with the surface of a DM dimensional 

1 /2 

sphere of radius Eq . This intersection will be also a 
sphere, which will be now of smaller radius and also of 
smaller dimension, D{M — 1). 

In order to compute ( [l4l| ), we change to the following 
notation 



V 
C 



((2Afi 



■ 1 Pm 1 Pi 

Nl/2 



y z 



■,Pm) 
(2MAf)^/^0,...,0,0,, 
(0,...,0,(2Afi)i/2,...,(2AfAf)'/',0,. 



,0) 
,0) 



= (0, . . . , 0, 0, . . . , 0, (2A/i)i/2, . . . , (2AfM)'/') (142) 



Note that these vectors satisfy Cx-Cy — Cy Cz ^ 
Cz-Cj: — 0. With these vectors so defined, Eqn. (141) 
becomes 



$(Po,So,Af) = J d''^'rs{r^ - E) d{Cx-v~PS) 



X 6{CyV-py) s{Cz-r-p^) 



D/2 



(143) 



Now we consider the following change of variables 

V ^r-{ -v^Pn + T^Pf! + T^f'n 1 (144) 



which is simply a translation and has unit Jacobian. Sim- 
ple algebra leads to 



M 



Uo 



X s{Cx-r') s{eyV') s{Cz-v') (145) 

where we have introduced the total internal energy 



Uo^iEo 



p2 

2Mo 



(146) 



where Aio = J2fj. ^^e total mass. 

We now consider a second change of variables 
A-P' through a rotation A such that 



AC, 


= (2Mo)'/'(l,. 


.,0,0,. 


.,0,0,. 


..,0) 


AC, 


= (2A^o)'/'(0,. 


.,0,1,. 


.,0,0,. 


..,0) 


AC, 


= (2A^o)'/'(0,. 


.,0,0,. 


.,0,1,. 


..,0) 



(147) 



It is always possible to find a matrix A that 
satisfies Eqns. (147). For example, consider 
a block diagonal matrix made of three identical 
blocks of size M x M. Then assume that each 
block is the same orthogonal matrix which trans- 
forms the vector ((2Afi)i/2^ (2A^2)^/^ ■ • • , (2AfM)^/^) 
into (2Alo)^^^(l, 0, . . . , 0). After the rotation (which has 
unit Jacobian and leaves the modulus of a vector invari- 
ant) the integral (|145D becomes 



M 



D/2 



Uo) 



X 5((2Mo)'/Vr) <5((2Xo)'/Vi'") S{{2Momn 



M „ 

= 1[{2M^)''^^ / rf^(^^-i)p"(5(7'"' - Uo) 



(148) 



We compute now the integral over the sphere in Eqn. 
( |148D by using that the integral of an arbitrary function 
F(x) — /(|x|) that depends on x only through its modu- 
lus |x| can be computed by changing to polar coordinates 



J F(x)rf^^x = UJM J f{r)r^''-^dr. 



(149) 



The numerical factor lum, which comes from the integra- 
tion of the angles, can be computed by considering the 
special case when f{r) is a Gaussian. The result is 
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UJM = 2 



rM/2 



r(M/2) 



(150) 



By using Eqn. (|l49| ), Eqn. ( |148| ) becomes 

M 

<i>(Po,i?o,M) = Y[i2AQ^/^ujDiM^i) 

X y dp S{p^ ~ U) (151) 

We need now the property 



(152) 



where are the zeros of /(x^) = 0. For the case of 
Eqn. (pll ) we have /(x) - [/, = ±({7)^/2 ^nd 

/'(x) = 2a:. Therefore, 



M 



$(Po,^o,M) = -cc;z5(Af-i)C/o 



(153) 



XIII. APPENDIX: VAN DER WAALS FLUID 

The free energy F of a system has as natural variables 
the number of particles N, the temperature T and the 
volume V, this is F = F{N,T,V). Its derivatives are 
given by the pressure, the entropy and the chemical po- 
tential 



P 
-S 
n = 



d£ 
dV 
dF 
df 
dF 
ON 



T,N 



V,N 



T.y 



(154) 



Because the free energy is a first order function of the 
extensive variables N,V , we have 



F{N,T, V)^Vf{n,T) 



(155) 



where n is the number densi ty an d f{n, T) is the fre e en - 
ergy density. The property (|l55| ) implies in Eqns. ( 154 ) 



-P{n,T) = f{n,T)-^ln 

« = -^KT) 

df , 
M=^KT) 

where s = S/V is the entropy density. 



(156) 



For a van der Waals gas, the free energy density is 
given by 

/(n,r) = -ksTn (^1 + ln ( ^ ) ) ) - (157) 

where the thermal wavelength is defined by 

h 



A(T) 



(27rmofcsT)i/2 



(158) 



Here, h is Planck's constant, mo is the mass of a molecule, 
and ks is Boltzmann constant. The constants a, b are 
the attraction parameter and the excluded volume, re- 
spectively. 

With Eqns. ( |l56| ) plus the well-known relationship / — 
e — Ts where e is the internal energy density, we easily 
arrive at the following relations 



e(n, s) — —kBT{n, s)n — an^ 



kBT{n, s) 



\2TTmoJ \1 



— nb 



2/D 



D + 2 



exp 



r 2s 

{DkBU D 



ksTn 2 

F(ri, s) = ; — an 

1 — no 

I nb f nb 

IJL[n, s) = ksT I r + In 



1 — nb \1 — nb 
2an 



A«(r) 



(159) 



which give the fundamental equation e(n, s) and the three 
equations of state in terms of the variables n, s. 

As it is customary for the van der Waals fluid, we in- 
troduce reduced variables 







f = 












p = 






a 




27b 


A = 




n = 


36n, 




36s 


s = 






ks 






e = 


^36e 




8a 



(160) 

In reduced units we have that the fundamental equa- 
tion becomes 



e(h, s) = — T(n, s)h n^ 

and the three equations of state are 

~, , 2rSD + 2 f n 

T[n, s) = cexp —<- h In 

D \ n 2 \6 — n 



(161) 
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P(n,S) = 

/ - \ 

(162) 

In this equation, we have introduced the the constant c 
which depends on microscopic parameters through 

27/,(^-2)/-D/,2 

c= = 4.836 x 10"^ (163) 

loTmioa 

where we have used the values corresponding to wa- 
ter in Z? = 3 dimensions, a ~ 1.5262 x 10^^*m^fcg/s^, 
6 = 5.0622 X 10-29m3, mo = 2.991 x IQ-^Sfeg, h = 
1.054 X 10"^'* Js Q. This constant appears in the the 
dimensionless quantity 



< - \ D/2 



(164) 



Note that in reduced units we still have the relation- 
ships 

^^~("'^~)^Mn,5) 



f(n,5) (165) 



dn 
de{n, s) 
ds 

Let us focus now on the distribution function (|llS| ) 
which in reduced units becomes 

P(n, S) = ^ exp{f |s - /3'='^*(?(n, S) - (166) 

where V = V jZh, /J'^^* = 8a^/27&, is_the reduced in- 
verse external temperature and — Xm^TIhl^a is the 
reduced external chemical potential. 

Due to the particular functional form of the funda- 
mental equation for the van der Waals gas, it is more 
convenient to study the distribution function of n, T in- 
stead of n, s. The Jacobian of the transformation is 
and, therefore, the new distribution function is 



P(n, T) = exp y{S(n, ~s) - ~s) - A^^'n)} 

(167) 
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